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Abstract 

This paper presents a model reduction method for the class of linear quantum stochas¬ 
tic systems often encountered in quantum optics and their related fields. The approach 
is proposed on the basis of an interpolatory projection ensuring that specific input-output 
responses of the original and the reduced-order systems are matched at multiple selected 
points (or frequencies). Importantly, the physical realizability property of the original quan¬ 
tum system imposed by the law of quantum mechanics is preserved under our tangential 
interpolatory projection. An error bound is established for the proposed model reduction 
method and an avenue to select interpolation points is proposed. A passivity preserving 
model reduction method is also presented. Examples of both active and passive systems are 
provided to illustrate the merits of our proposed approach. 


1 Introduction 

Over the past few decades, considerable interest has been drawn to quantum systems involv¬ 
ing open harmonic oscillators linearly coupled to one another and to external Gaussian fields, 
especially in areas such as quantum optics, optomechanics, and superconducting circuits. This 
type of systems is used to model quantum devices such as optical cavities, mesoscopic mechan¬ 
ical resonators, optical and superconducting parametric amplifiers, and gradient echo quantum 
memories (GEM); e.g.. see [Tm. A number of applications of linear quantum stochastic sys¬ 
tems have been theoretically proposed or experimentally demonstrated in the literature. In 
particular, they can serve as coherent feedback controllers [^|^, i.e., feedback controllers that 
are themselves quantum systems. In this context, they have been shown to be theoretically ef¬ 
fective for cooling of an optomechanical resonator , can modify the characteristics of squeezed 
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light produced experimentally by an optical parametric oscillator (OPO) [^, and, in the setting 
of microwave superconducting circuits, a linear quantum stochastic system in the form of a 
Josephson parametric ampliher (JPA) operated in the linear regime has been experimentally 
demonstrated to be able to rapidly reshape the dynamics of a superconducting electromechani¬ 
cal circuit (EMC) [^. Linear quantum stochastic systems can also be used as optical Liters for 
various input signals, including non-Gaussian input signals like single photon and multi-photon 
states. As Liters they can be used to modify the wavepacket shape of single and multi-photon 
sources [^[^. Also, linear quantum stochastic systems can dissipatively generate Gaussian 
cluster states [10] as an important component of continuous-variable one way quantum comput- 
11 . They can also be exploited for classical signal processing on quantum devices, such as 


ers 


in light processors where photons transport information between cores on a chip and between 
chips [I2] . 

The linearly coupled open harmonic oscillators can be completely represented in Heisenberg 
picture of quantum mechanics by a class of linear quantum stochastic systems described using 
quantum stochastic differential equations (QSDEs) and subsequently, a quartet of matrices 
{A, B,C, D) (analogous to classical non-quantum linear systems) [5,13-17. However, unlike 
many classical non-quantum systems, these matrices cannot be arbitrary and they are not 
independent of one another due to the restrictions imposed by the laws of quantum mechanics. 
In fact. A, B, and C must satisfy certain equality constraints while D must be of a speciLc 
form. These constraints are referred to as the physical realizability conditions [^. The class of 
linear quantum stochastic systems provides us with a tool to develop fundamental principles 
in quantum control theory in a similar way that classical non-quantum linear control systems 
have facilitated the advancement of classical systems and control theory. 

In many linear quantum control problems, the input-output relation of the system (or con¬ 
troller) described by its transfer function is more important than its realization described by 
the matrices {A, B,C, D) e.g., see [^[6[ [I^[l8[|^ . Unfortunately, the system transfer functions 
may be complex involving a large degree of freedom and hence, physically implementing systems 
with such transfer functions can be a challenging task. In these situations, it is more appealing 
to construct an approximating system with a lower degree of freedom such that the transfer 
functions of the original and the approximating systems are closely matched. The process of 
constructing a lower order approximating system is known as model reduction. The contribution 
of this paper lies in this model reduction problem. 

A challenge facing model reduction in quantum setting is the retention of the physical 


realizability property. In |^ -25 , singular perturbation methods (also known as adiabatic elim¬ 
ination) were studied for reduction of linear quantum systems comprising subsystems that 


evolve at well-separated time-scales. In 26 , an eigenvalue truncation method was proposed for 
a sub-class of completely passive systems with distinct poles. More recently, an adaption of 
the well-known balanced truncation method was proposed in [27] for linear quantum stochastic 
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systems whose controllability and observability Gramians satisfy some commutation conditions. 

As the existing model reduction methods are limited to sub-classes of systems with specific 
properties, this paper presents a model reduction approach that allows an approximation of a 
more general linear quantum stochastic system. In particular, we propose a model reduction 
method on the basis of the tangential interpolatory projection method introduced in [28] . This 
approach constructs an approximating system such that its transfer function interpolates the 
transfer function of the original system at multiple points along some tangent directions. Impor¬ 
tantly, we show that the reduced system is physically realizable. An error bound is established 
for the proposed method, and an avenue to select interpolation points and tangent directions 
are presented. 

The remainder of this paper is structured as follows. In Section we define the class of 
linear quantum stochastic systems under consideration. The tangential interpolatory projection 
model reduction approach is proposed in Section Section presents a passivity preserving 
model reduction method. Concluding remarks are then provided in Section 


2 Preliminaries 


2.1 Notation 


We will use the following notation: i = * denotes the adjoint of a linear operator as well 

as the conjugate of a complex number. If A = [ajk] then and = (A^)"*", where 

(•)"'' denotes matrix transposition. 3f?{A} = {A + A^)/2 and ^{A} = ^(A — A^). We denote 
the identity matrix by I whenever its size can be inferred from context and use In to denote an 
nxn identity matrix. Similarly, Omxn denotes a m x n matrix with zero entries but we drop the 
subscript when its dimension can be determined from context. We use diag(Mi, M 2 ,..., M„) 
to denote a block diagonal matrix with square matrices Mi, M 2 ,..., Mn on its diagonal, and 
diag„(M) denotes a block diagonal matrix with the square matrix M appearing on its diagonal 


blocks n times. Also, we let J = 


0 1 

-1 0 


and In 


In I = diag„(JI). We use N to denote 


the set of natural numbers. We use e, = [0,0,..., 0,1,0,..., 0]"'" to denote an indicator vector 
with one in the i-th position and zero elsewhere. For a subspace H of a vector space, we write 
P-H to denote an orthogonal projection onto P. We also write to denote the orthogonal 
complement of P. 


2.2 The class of linear quantum stochastic systems in quadrature form 

A linear quantum stochastic system [5T3, 17 is an open Markov quantum system involving n 
independent quantum harmonic oscillators coupled to m independent external continuous-mode 
bosonic fields. Let x = {qi,pi,q 2 ,P 2 , ■ ■ ■, Qn,Pn)~^ denote a vector of the canonical position and 
momentum operators of a n degree of freedom quantum harmonic oscillator satisfying the canon- 
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ical commutation relations (CCR) xx^ — {xx~^)^ = 2iSn- That is, qj and pj are position and 
momentum operators of j-th harmonic oscillator, respectively. Let Y{t) = ... ,Ym{t))~^ 

denote a vector of continuous-mode bosonic output fields that results from the interaction of the 
quantum harmonic oscillators and the incoming continuous-mode bosonic quantum fields in the 
m-dimensional vector A{t). Note that the dynamics of x{t) is linear, and Y(t) depends linearly 
on x(s) with 0 < s < t. Following [^, the dynamics of a linear quantum stochastic system can 
be described in quadrature form as 


dx{t) = Ax{t)dt + Bdw{t); x(0) = x. 

dy{t) = Cx{t)dt + Ddw{t), (1) 

where A € M 2 nx 2 n^ ^ ^ ^ 2 nx 2 m^ Q ^ ^ 2 mx 2 n^ ^ ^ 2 mx 2 m^ 

w{t) = 2(3f?{^i(t)}, 9={^i(t)}, 3={^2(t)}, • • • , 

y{t) = 2[^{Yi{t)}Myi{t)}MY2{t)}My2{t)},---Myn.{t)}Myrn{t)}y. 


Here, w{t) (input fields in quadrature form) is taken to be in a vacuum state where it satisfies 
the Ito relationship dw{t)dw{t)~^ = {I + tSm)dt] see [^. Note that in this form it follows that 
is a real unitary symplectic matrix. That is, it is both unitary (i.e., DD^ = D = I) and 
symplectic (i.e., D^rnD^ = Jm)- However, in the most general case, D can be generalized to 
a symplectic matrix that represents a quantum network that includes ideal squeezing devices 


acting on the incoming field w{t) before interacting with the system 15,17. In general one may 
not be interested in all outputs of the system but only in a subset of them, see, e.g., [^. That 
is, one is often only interested in certain pairs of the output field quadratures in y{t). Thus, in 
the most general scenario, y{t) can have an even dimension 2i < 2m. 

Unlike classical non-quantum systems, the matrices A, B, (7, H of a linear quantum stochas¬ 
tic system cannot be arbitrary and are not independent of one another. In fact, for the system 
to be physically realizable [5||17|, meaning it represents a meaningful physical system, they must 
satisfy the constraints (see [^[I^[l7l[^ ): 


AJji ^nA^ = 0, 

(2) 

-|- = 0, 

(3) 


(4) 


3 Tangential Interpolatory Projection for Model Reduction 


Interpolatory projection model reduction framework was developed by De Villemagne and Skel¬ 


ton in 29 to construct a reduced-order system such that its transfer function interpolates the 


transfer function of the original system at a single point. Limitations of this single-interpolation- 
point approach were discussed in pOpT] including the loss of approximation accuracy away from 
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the single interpolation point. For this reason, multi-point interpolatory projection framework 
was introduced. In particular, two key tangential interpolatory projection approaches were pro¬ 
posed in 28 : the left and the right tangential interpolatory projections. Let H(s) and Er{s) 


be the transfer functions of the original and the reduced-order systems, respectively. Given 
a set of 2r interpolation points {cri,a 2 , ■ ■ ■ ,cr 2 r} and 2r left (or output) tangent directions 
• ■ ■, lJ- 2 r} C the left tangential interpolatory projection approach ensures that 


nl'^ricTi) = 


(5) 


for all i = 1, 2,..., 2r. Similarly, the right interpolatory approach ensures that 


= E{ai)ui ( 6 ) 

for all i = 1, 2,..., 2r, where i^i, 122 ,, V 2 r £ are the right (or input) tangent directions. 

We stress that we are interested in constructing a reduced system while meeting either ([^ or 
Q. In general, one might consider the left tangential interpolatory projection approach (instead 
of the right) because the left tangent directions belong to a smaller complex space than (or 
at most equal to) the space where each right tangent direction lives in (since i < m for the 
class of linear quantum stochastic systems). The smaller space reduces the choices of tangent 
directions and may simplify the selection of appropriate tangent directions. Intuitively, this may 
also help in mitigating the effects the use of inadequate tangent directions may have on model 
approximation error. However, in some situations, it may be more appropriate to consider the 
right tangential interpolatory projection e.g., when one wants to place more emphasis or weights 
on the input components that are of interest. For this reason, we will present both the left and 
the right tangential interpolatory projection approaches. 


3.1 Petrov-Galerkin Projection Approximation 

The tangential interpolatory projection method can be achieved by constructing a reduced-order 
model of a linear quantum stochastic system via Petrov-Galerkin projection approximation. 
That is, given a full-order linear quantum stochastic system in quadrature form ([^, we construct 
a reduced-order linear quantum stochastic system of the form 

dXr{t) = ArXr{t)dt + Brdw{t); Xr(0) = X 

dyr{t) = CrXr{t)dt + Drdw{t) (7) 

where Ar = wj AVg, Br = wj B, Cr = CVg, Dr = D. Here, Wq G M2nx2r ^ ^2nx2r 

the left and right projection matrices, respectively, which satisfy the condition Vq = I. 

Let H(s) = C {si — A)~^ B and Er{s) = Cr {si — Ar)~^ Br. Since we consider Dr = D, 
it can be seen that the left tangential interpolation condition (|^ can be achieved by meeting 
a simpler condition /it“((jj) = /j,\3r{ai) for all i = l,2,...,2r. Similarly, the right tangential 
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interpolation condition Q is attained when 'E{ai)vi = is satisfied for alH = 1,2,..., 2r. 

On the basis of these simplified interpolation conditions, we can now re-state an important 
result on tangential interpolatory projection in terms of our linear quantum stochastic systems 
in quadrature form and Petrov-Galerkin projection approximation. 

Theorem 1. l3^ Theorem 1] Given a full-order linear quantum stochastic system written in 
quadrature form 0 with the transfer function H(s), consider a reduced-order model of the form 
Q constructed through Petrov-Galerkin projection approximation with the transfer function 
Er{s). Then we have that 

(i) the left interpolation condition ([^ holds if Wq has full rank and £ 

range(VFq) for all i = 1,2,, 2r. 

(ii) the right interpolation condition Q holds if Vq has full rank and {ail — £ 

range(V"g) for all i = 1,2,..., 2r. 

For many classical (non-quantum) systems, it is possible to choose Wq and Vq independently 
so that both the left and right interpolation conditions 0-0 are satisfied. Unfortunately for 
quantum systems, due to the law of quantum mechanics, the projection matrices Wq and Vq 
cannot be independently selected. Hence, using the above result, we will present avenues to 
choose Wq and Vq such that either the left interpolation condition 0 or the right interpolation 
condition 0 is satisfied, while ensuring that the physical realizability property of a linear 
quantum stochastic system is preserved. However, before presenting our physical realizability 
preserving model reduction method, let us recall the following well-known result. 

Lemma 2. For any full-rank (non-singular) real skew-symmetric matrix 0 £ (i.e., 

0 = —0"''^, there exists a non-singular (full-rank) real matrix T £ such that TQT~^ = 

Jn. 

3.2 Proposed Physical Realizability Preserving Model Reduction 

Physical realizability is an important property of a quantum system. In constructing a reduced- 
order model of a linear quantum stochastic system, it is important to ensure that the reduced 
system represents a meaningful physical quantum system. As previously mentioned, the phys¬ 
ically realizability property places some restrictions on the projection matrices Wq and Vq. 
Motivated by the method presented in j^, we will now present avenues to choose Wq and Vq 
so that the reduced system is guaranteed to be physically realizable. 

3.2.1 Left Tangential Interpolatory Projection 

To achieve the left tangential interpolation condition 0 , consider a subspace 

Wq = span!. (plC {ail - A)~^y\ (8) 

I ^ ^ J i=l,2,...,2r 
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where interpolation points, {cii, cJ 2 , ..., cr 2 r}, and the left tangent directions, {fii, fj, 2 , ■ ■ ■, f^ 2 r}, 
are chosen such that 

(i) (ail — A) is invertible for alH = 1, 2,..., 2r, 

(ii) / 0 for alH = 1, 2,..., 2r, 

(hi) the subspace Wq has dimension 2r, and 

(iv) a real basis Wq of Wq exists (see Remarkon the existence of a real basis) and WjSnWq 
has full rank. 


We highlight that the subspace is required to have dimension 2r to ensure that the the reduced 
system has r degrees of freedom. We also note that WjInWq has full rank if and only if the test 

0 .; 


matrix 


WJ 


InWa 


has full rank; see 


the test matrix would generally have fuf 


34 


Eq. (3.4)]. From direct inspection, we can see that 


rank except for some Wq with specific structures. In 


other words, {cri, (T 2 , • ■ •, o' 2 r} and {m, ^ 2 , ■ ■ ■, ^J' 2 r} can typically be chosen such that W^ InWq 
has full rank except for systems with specific forms of A and C. 

We then propose the projection matrices, Wq and Vq, to be 


Wq = WqT^ (9) 

Rq = SnWq (Wj JnW?) (10) 

where T G ]^ 2 rx 2 r jg non-singular transformation matrix such that T{WjInWq)T~^ = The 
existence of T is guaranteed by Lemma when Condition (iv) stated above holds. Here, Wq 
is also a real basis of Wq because T is a non-singular matrix. Thus, the product WjSnWq is 
invertible because WjSnWq is assumed to have full rank. We highlight that, by choosing the 
above projection matrices, Wj Vq = I as required by Petrov-Galerkin projection approximation. 
We now present a new result showing that, using the above projection matrices, our proposed 
reduced-order linear quantum stochastic system Q is physically realizable while meeting the 
interpolation condition ([^. 


Theorem 3. Given a physically realizable full-order linear quantum stochastic system written 
in quadrature form 0 with the transfer function E(s), consider a reduced-order model of the 
form Q with the projection matrices Wq and Vq given by 0 and respectively. Then 

the reduced-order system is also physically realizable and its transfer function Erfs) interpolates 
H(s) in the sense of 0. 

Proof. First note that, because Wq is a basis of the subspace Wq defined in ([^, Wq has full 

rank and _ 

^/xtC(cJi/ — H)“^^ G range(ITq). By Theorem 


Sr(s) interpolates S(s) in the sense of ([^. 
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Now we will show that the reduced-order system with the transfer function is phys¬ 

ically realizable. Let us first recall that Wq = WqT~^ where T is an 2r x 2r matrix such that 


T{Wg SnWq)T~^ = Sr- Hence, we have that W'SnWq = Sr- It then follows from (10) that 


SnWq = VqSr and Sn = SrV^ . Using these identities, pre- and post-multiplying the first 
physical realizability condition of the full-order system ([^ by and Wq, respectively, gives 


us 


Wj ASr^Wq + Wj SnA^Wq + BSmEA Wq = 0 


ArSi 


-|- BrSmB^ — 0 . 


Similarly, pre-multiplying the second physical realizability condition of the full-order system 
© by Wj gives us that SrCj + BrSmD~^ = 0. The theorem statement then follows because 
Dr = D. m 


From the above theorem, we stress that the proposed projection matrices ([^-([Tol) lead to 
a physically realizable approximate system. Our proposed tangential interpolatory projection 
method can be applied to a large class of linear quantum stochastic systems compared to existing 
methods proposed in |27[|35| , as will be demonstrated in examples. Moreover, the tangential 
interpolatory projection framework provides an appropriate reduction scheme when we are only 
interested in specific input-output responses of the system at a particular range of frequencies, 
Recall, for a subspace H of a vector space, that we write to denote an orthogonal 
projection onto Ti. Let X C and y C We use (p iX,y} to denote the largest principal 
angle between the subspaces X and y, which is defined as 

cos{(j){X,y)) = inf \{x,y)\- (11) 

xGX ,y^y 

ll*ll=ij|y||=i 


Note that cos {4> {X,y)) = -^1 — \\Px — Py P 36 . We now present an error bound for the left 
tangential interpolatory projection method. 


Proposition 4. Consider a reduced-order model constructed using the subspace Wg defined in 
(|^, and the full-rank projection matrices Wq and Vq as given in (© and ( [Tol ), respectively. Let 
hlw{s) = ker{Wg’'(s/ — ^)} and Uv{s) = ker{I7^(s/ — ^)I}. For any s G C that is not an 
eigenvalue of either A or Ar, we have that 


||H(s) - Er{s)\\ = \\C{sI - A)-^ (/ - Q(s)) P|| , (12) 

= ||C(/-P(s))(s/-A)-ip||, (13) 

where Q{s) = {si — A)Vq{sI — Ar)~^Wq and R{s) = Vq{sl — Ar)~^Wq {si — A) are oblique 
projection operators (i.e., Q{s) and R{s) are idempotent), and when A and Ar have no poles 
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on the right half plane and the imaginary axis, we have the following %oo bounds: 


■"J’lloo 


< sup ( (l - ||Pw± - Puyi^u,)\^ 
X C{lUjI — A)~^Py^;± 

< sup ( (l - \\P-v± - ) 

ujGR V ^ ^ ^ 


- 1/2 


PUvil^^)^! 


- 1/2 


X CP; 


Viwi'^^) I 


P^riwol - A)-^B\\y 


(14) 


(15) 


Proof. Let us first recall that the range and kernel of a projection operator P are complementary 
in the sense that ker{/ — P} = rangejP} [^. Under our assumption on s G C, ker{(2('S)} = 
range{/— (5(s)} = because {si — A)Vq{sI — Ar)~^ is full rank. Note that ker{/ —(5(s)}-*- = 
range(/—Q(s)^^)} |37l Theorem 6.6]. From this identity, the fact that Wj{sI—Ar)~^ is full rank, 


and the definition of (10), we have that ker{/— Q(s)}-*- = ker{(5(s)'l'} = ker{Vg^ (si — A)l} = 
Uv{s). 

From the above identities, we have that I — Q{s) = Pvyv {I — Q('S)) Pu{s)- Thus, we have 
that 


“(s) - “^(s)|| 


C{sl-A)-^ 


I - (si - A)Vq{sI - Ar)-^Wj 


B 


C{sl - A)-^ {I - Q{s)) B\\ 

C{sl - A)-^P^r {I - Q{s)) Puyis)B 


This establishes 


. The result (14) then follows from the definition of the Poo norm and that 

- 1/2 


36 


Theo- 


||/-Q(s)|| = sec(range{/-Q(s)},ker{/-Q(s)}-L) = {l - \\Pwf - Puiiuj)^) 
rem 6.1]. Finally, ( |13[ ) and ( |15[ ) follow from similar arguments to the above. This establishes 
the proposition statement. ■ 


Note that for any s G C, Pw^j Pvf-: Puw{s)^ Puv(s) can be computed from 
{/ii, /i 2 , • ■ ■, /i 2 r}) {ui, (T 2 , ..., o' 2 r}, Jn, A, and C. Thus, the bound (14) can be obtained without 
computing Wg, Vg, and the reduced system model. 

Moreover, we note that the two Poo bounds, (14) and ( [I^ , are established on the basis of 
the operator norms of two different oblique projection operators Q{s) and R{s). Since they are 
oblique projection operators, their operator norms may be different and large (unlike orthogonal 
projections). Thus, the two upper bounds can be different and one might be tighter than the 
other, depending on {/ii,// 2 , • • ■ ,h 2 r}, Wi,o' 2 , ■ ■ ■ ,cr 2 r}, Jn, A, and C. 


Remark 5. As we require Wg and Vg to he real-valued matrices, the interpolation points cannot 
he any arbitrary complex numbers and the tangent directions cannot he any arbitrary complex 
vectors. In fact, a real basis of the subspace Wg exists when the interpolation points, ai,..., a 2 r, 
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and their corresponding tangent directions, /Ui,... ,ti 2 r, o,re real or occur in conjugate pairs; see 


also 38. Remark 


41- 


3.2.2 Right Tangential Interpolatory Projection 

The right tangential interpolatory projection follows similar construction to the left tangential 
interpolatory projection presented previously. It involves a different subspace which is defined 
as 


Vq = span{(cri/ - A) ^Bvi] 

where interpolation points ai,a 2 , ■ ■ ■ ,(T 2 r and the right tangent directions vi,n 2 , ■ ■ ■ ,n 2 r are 
chosen such that: (i) {ail — A) is invertible for all i = 1,2, ...,2r, (ii) t'j 7 ^ 0 for all i = 
1 , 2 ,..., 2 r, (hi) the subspace Vq has the dimension of 2 r, and (iv) a real basis Vq of Vq exists 
and Vq ^nVq has full rank. We then propose Vq and Wq to be 

Pg = %T'^ (16) 

Wq=^nVq{yJ^nVqY^ (17) 

where T G ]^ 2 rx 2 r jg non-singular transformation matrix such that and 

Vq is a real basis of Vq. 


Theorem 6. Given a physically realizable full-order linear quantum stochastic system written 
in quadrature form Q with the transfer function H(s), consider a reduced-order model of the 
form Q with the projection matrices Wq and Vq given by (16) and respectively. Then 

the reduced-order system is also physically realizable and its transfer function Erfs) interpolates 
H(s) in the sense of (©. 

Proof. This proof follows similar argument to the proof of Theorem [fusing the result of Theorem 

0ii)- ■ 


3.3 Selection of interpolation points and tangent directions 

In this subsection, we propose a heuristic to select interpolation points and left tangent directions 
for our reduced-order linear quantum stochastic system. Note that we will only consider the left 
tangential interpolatory projection approach. However, similar approaches can be undertaken 
for the right tangential interpolatory projection approach. 

For the left tangential interpolatory projection, let us assume that the output fields y{t) 
can be re-ordered from the most important pair (of momentum and position operators) to 
the least through a permutation matrix Hy. We now present an approach to heuristically 
choose left tangent directions so that the more important output fields are matched at a larger 
(or at least equal) number of frequencies. Recall that Ci = [0,0,..., 0,1, 0,..., 0]"'" G 
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denotes an indicator vector with one in the i-th position and zero elsewhere. Also let Ef^ = 
(ei, ei, 62 , 62 ,, Cfc, Cfc) for any k £ N. We propose the tangent directions to be 

(/il,^2, • •. ,M2r) 


= nj X 


Er, if r < £, 

{Ee,..El, 61 , 61 , 62 , 62 , ■ ■ ■), otherwise. 


(18) 


Here, the tangent directions are chosen in conjugate pairs to ensure that a real basis of Wg 
exists. 

It now remains to choose the interpolation points. In this approach, we are interested in 
matching frequency responses of the full-order system. Therefore, we propose that the set of 
interpolation points be along the imaginary axis of the form 


((Ti,o- 2, ... ,a2r) = (^a;^,-ZW2, • ■ • (19) 

where wj,... > 0. Again, the interpolation points are chosen in conjugate pairs to ensure 

that a real basis of Wq exists. For the chosen tangent directions and any H = (uji,uj 2 , ■ ■ ■, ^r) £ 
M'’ (such that the conditions under ([^ hold), let VF(H) be an orthonormal basis of the subspace 

span |(za;j/-Al)“^CVfc, 

where k = 2j - 1. Let H(H) = J„IT(H)(IT(H)tj„IT(H))-i and 
Q{s, n) = {si - A)v{n){si - w{nyAV{n))-^w{ny. 


3.3.1 Hoo based selection of interpolation points 


Since we are interested in matching frequency responses of the full-order system, it is natural to 
attempt to minimize the error between H(s) — Hj,(s) when s is purely imaginary (i.e., minimize 
error across all frequencies). Thus, we define a cost function on the basis of the Hoo norm based 
on the error formulae © as 


Joo(fl) 


sup 

ojeK 


C{iujl - A)-^ 


(^I - Q{iuj,n)^ B 


The vector of interpolation points = {oj 1 ,U 2 , ■ ■ ■ ,u)^) is then chosen such that it is a local 
minimizer of 


3.3.2 ?^2 based selection of interpolation points 

For some quantum systems, it is possible that the cost i7oo(^) is the same for all H G M'’; 


see e.g.. Example 13 In such cases, it may be more appropriate to consider an optimization 
problem based on ?^2 error. Let E{uj, Q) = C{iujl — A)~^ — Q{vjJ, H)^ B Let us define a cost 

function on the basis of the I-L 2 norm based on the error formulae (12) as 


J 2 { 11 ) = J trace (^E{u,Q,)^E{uj,Q)'j dio. 
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The vector of interpolation points ■ ■ ■ ,u}^) is then chosen such that it is a local 

minimizer of c 72 (f^)- 


3.4 Illustrative examples (active systems) 

In this subsection, we present some practical examples illustrating the application of our pro¬ 
posed tangential interpolatory projection approach on active linear quantum stochastic systems 
(having some squeezing). Particularly, the systems considered in these examples cannot be trun¬ 
cated through the existing quasi-balanced truncation method |27| because the controllability 
and observability Gramians of the full-order systems do not satisfy the commutation condition 
for the existence of a symplectic balancing transformation matrix. 


Example 7. Consider an optomechanical system proposed in 39l for back-action evading mea¬ 


surement of position, consisting of an optical cavity with two movable mirrors, conceptually il¬ 
lustrated in Tig. For an introduction to optomechanical systems we refer the reader to 40,4^, 
while for studies of control on such systems see, e.g., The cavity is pumped by a strong 

coherent laser and each mirror is subjected to radiation pressure and thermal noise. This system 
has three degrees of freedom comprising an oscillator inside the cavity, described by the quadra¬ 
tures {qi,pi), and two mechanical oscillators from the motion of the two mirrors described by the 
quadratures {q 2 ,P 2 , Q 3 ,P 3 )- The dynamics of this optomechanical system can be linearized about 
the steady-state value of the mean amplitude of the cavity mode, see, e.g., for details. The 

linear approximation can be described in the quadrature form 0 with X = {qi,Pi,q 2 ,P 2 , 93,93^ 


and the following system matrices 391: 


A = 


_ K 

2 

0 

0 

-r 

0 

0 


0 

0 

0 

0 


0 

-r 

_7 

2 

0 

0 

—Q/, 


0 

0 

0 

_2 

2 

Qh 

0 


0 

0 

0 

—rifc 

_7 

2 


0 

0 

Qh 

0 

0 


0 


B = diag{^/KI 2 , C = [-v/k /2 02 x 4 ], and D = [—I 2 02 x 4 ]? where k > Q is the cavity 

decay rate, 'y > 0 is the damping rate of the two mechanical oscillators, P > 0 is the optome¬ 
chanical coupling rate (due the coupling between the mirror degrees of freedom and the cavity 
mode via radiations pressure), and idb is the system half-bandwidth. Here, inputs 1-2 are the 
two quadratures of the laser field while inputs 3-6 describe the thermal fluctuations acting on 
the mirrors. Note that this system is active as the coupling P leads to a squeezing Hamiltonian. 

Let K = 2 X 10^, 7 = 100, P = 7.0711 x 10^, and LI = 10^. The poles of the system are 
then at —50 ± lO^z and —10®. This choice of parameter values allows us to compare our model 


reduction method with the singular perturbation approximation presented in (25 1. 
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Figure 1: 



Conceptual model of the optomechanical system. 



Figure 2: Example 0 The frequency responses (Bode plot) of the full-order linear quantum stochastic 
system (black solid line) and the tangential interpolatory projection approximation (red dashed line). 
The center frequency (where the peaks appear) is the average resonant frequencies of the two mechanical 
oscillators, i.e., {ujmi +Wm 2 )/ 2 . 


Consider an approximation with r = 2. We are interested in matching frequency responses 
from the thermal fluctuations on the mechanical oscillators to the quadratures of the output 
field Yi{t). We apply the right tangential interpolatory projection model reduction approach to 


provide different weights on the inputs. As proposed in (18), we choose the tangent directions to 


be 1 ^ 2 , 1 '^, Vi) = ( 65 , 65 , 66 , 66 ). The corresponding interpolation points of the form (19) are 
obtained by finding a local minimizer of the TLoo based optimization problem. For computational 
simplicity, we set cuf = With this assumption, we find = 1.05 x 10^ as a local 

minimizer. We note that this leads to a stable reduced system whose poles are at —50 ± 10^^. 

Fig. d illustrates the frequency responses of the original and our reduced systems at output 
2, using black solid and red dashed lines, respectively. Note that the responses from inputs 1, 
4, and 5 are omitted as the transfer functions from these inputs are zero due to the dispersive 
coupling. Similarly, the responses at output 1 are also omitted because, for both systems, the 
transfer functions from inputs 2-6 are zero and the ones from input 1 are the same as those at 
output 2 from input 2 (shown in the figure). From the figure, we see that the responses of the 
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Figure 3: Example The approximation error ||(.:i(s) — ,:ir(s))||- 



Figure 4: Examplej^ The approximation error at each purely imaginary point (or frequency) ||(^(^a;) — 
Sr(*w))||. 


two systems are well matched over the narrow bandwidths (2Qb) of the forces (inputs 3 and 6) 
acting on the mechanical oscillators. The magnitude responses from the pump beam (input 2) of 
the two systems are the same, whilst the phase response of our reduced system matches that of 
the original system at high frequencies. The error ||(H(s) — Hr(s))|| for each s is shown in Fig. 
dandg} Note that the tangential error is large when s is close to the poles of the systems. The 
TLoo error of the reduced system is 2.00. The TLoo bounds (14) and (15) are 2.45 and 3.96 x 10^, 


respectively. Thus, we see that the bound (14) is tight, while the bound (15) is too conservative 
for this example. 


For comparison, singular perturbation approximation 25l is applied to eliminate the cavity 
oscillator {qi,pi) and the TLoo error of the singular perturbation approximation is 4.45. Flence, 
our model reduction method provides a better reduced model than the singular perturbation ap- 
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proximation (in terms of T-Loo error) for this optomechanical example with the given set of 
parameter values. 

Example 8. Consider a quantum control system originally considered in Section IV], com¬ 
prising of a cascade connection of an optical parametric oscillator (OPO), two optical cavities, 
and a stabilizing linear quantum controller. Let Uq{t) and Up{t) denote the momentum and po¬ 
sition quadratures of the output of the stabilizing controller, respectively. The dynamics of the 
overall control system is described by 


dx{t) = Ax{f)dt + Budu{t) + Bdw{t) 
dy{t) = Cx{t)dt + Ddw{t) 


where du{t) = [uq{t),Up{t)]^dt + dv{f) for some vacuum quantum noise vector v{t) different 
from w{t) originating from the controller output field, and 


A = 


-0.5006/2 -0.0374/2 -0.0410/2 

0 A22 -1.0954/2 , 

0 0 -0.6/2 


A22 — 


-1 -1.05 

-1.05 -1 


B,, = 


n T 


- 0 . 0374/2 -I2 - 1 . 0954/2 


c = 


I 2 


02x4 



D = 


I 2 



Here, the plant is unstable (the eigenvalues of A are 0.05, —0.05, —0.5007, and —2.05). 

To stabilize the plant, we first design a classical (non-quantum) output feedback controller 
using the pole placement method. We select the closed-loop poles to he at —0.2, —0.3, —0.5, 
—0.6, —0.9, —1.5, —1.8. The controller dynamics is then described by 


dz{t) = Acz{t)dt + Bcdy{t) 
du{t) = Ccz{t)dt 
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where 


■ - 1.5500 

- 0.0001 - 

- 0.0016 

0.0000 - 

0.0139 

o . oooo ' 

- 0.0052 

- 2.4500 

0.0000 - 

- 0.0315 

0.0000 - 

0.0447 

- 10.3270 

- 0.0022 - 

- 1.0920 - 

- 0.0002 - 

0.3718 

0.0004 

0.2787 

39.3723 

0.0003 

0.2090 

0.0001 - 

1.1943 

17.5312 

0.0014 

1.0494 - 

- 0.0002 

0.1927 

0.0005 

_ - 0.0207 

0.0007 

0.0003 

0.1741 

0.0002 - 

0 . 7083 _ 

1.0493 

o.ooof 


■- 0.0006 

o . oooo ' 

T 

0.0052 

1.9493 


- 0.0000 

- 0.0007 


10.3276 

0.0022 

, Cc = 

- 0.9580 

- 0.0003 


- 0.2787 

- 39.3717 


0.0002 

- 0.1590 


- 17.5305 

- 0.0014 


- 0.7236 

- 0.0001 


0.0207 

0.0000 


- 0.0004 

0.0989 



This type of controller can be made physically realizable as a linear quantum system with six 
degrees of freedom by introducing seven additional inputs v{t) to the controller ^ Lemma 5.6]. 
That is, there exists a matrix By such that the system described by 

dz{t) = Acz{t)dt + [By B^[dv{ty' dy{t)~^]^ 
du{t) = Ccz{t)dt + [I 2 0][d?;(t)''~ dy{t)~^]^ 

is physically realizable (as a quantum system). Using the formula given in ^ Lemma 5.6], 
By G is given by 


0.0007 

0.0000 

0.1590 

- 0.0003 

- 0.0989 

- 0.0001 

- 0.0000 

0.0006 

0.0002 

0.9580 

- 0.0004 

0.7236 

- 0.0328 

0.0000 

0.1957 

- 0.0012 

- 0.3337 

0.0004 

- 0.0000 

- 29.6894 

- 0.0002 

- 0.1134 

0.0003 

0.0001 

0 

- 0.0328 

- 0.0000 

- 0.0387 

0.0000 

- 0.0002 

29.6921 

0.0000 

- 0.0056 

0.0000 

0.0096 

- 0.0000 

- 0.0387 

0.0012 

- 8.1505 

- 0.0000 

13.8054 

- 0.0163 

- 0.0000 

- 0.1134 

- 0.0001 

- 24.9918 

- 0.0001 

- 0.0027 

0.0000 

0.1957 

- 0.0000 

- 8.1505 

0.0005 

- 0.0006 

- 0.0056 

0.0002 

28.4766 

0.0001 

2.0589 

- 0.0024 

- 0.0002 

- 0.0004 

- 0.0006 

0.0163 

- 0.0112 

0 

0.0000 

0.0001 

0.0024 

- 0.0027 

- 0.0041 

- 29.6921 

- 0.0000 

- 0.3337 

- 0.0005 

13.8054 

0.0000 

- 0.0112 

0.0096 

- 0.0003 

2.0589 

0.0001 

26.2046 

0.0041 


Note that By does not affect the location of the poles of the closed-loop system; see ^ Eq. 
( 24 )]. Also, note that the controller is stable with the poles at —0.3837, —0.4597 ± —0.6664^, 
-0.7930, and -1.2726, -2.0299. 

In this example, we are interested in approximating the controller described by 
{Ac, [By Be], Cc, [I 2 0]) with r = 2. To place more emphasis on the input dy{t) of the controller 
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Figure 5: Example The frequency responses (Bode plot) of the full-order controller (black solid line) 
and the tangential interpolatory projection approximation (red dashed line) from inputs 15 and 16 (i.e., 
output of the plant y{t)). 


(output of the plant), we consider the right tangential interpolatory projection method. We 
choose (z^i, 1 ^ 2 , 2 ^ 4 ) = (ei 5 , ei 5 , ei 6 , eie). The corresponding interpolation points of the form 

(19) are obtained from a local minimizer of the Tioo optimization problem described in the pre¬ 
vious subsection. Again, we simplify the problem by assuming that = UJ 2 = oj‘^. With this 
assumption, we have that = 0.2900. This leads to a stable reduced controller with the 
poles at -0.2576 ± 1.4795Z, -0.5391, and -1.4958. 

Fig. 0 illustrates the frequency responses of the original controller and the tangential in¬ 
terpolatory projection approximation from the input y{t). We see that the responses at output 
1 from input 15 is closely matched across all frequencies. The other pairs of input-output re¬ 
sponses from the quadratures in y{t) to quadratures of the controller output u{t) are also quite 
well matched except around the frequency of 0-2 rad/s. The closed-loop system with the approx¬ 
imate controller has poles at—0.1265±0.1404z, —0.3272, —0.3654± 1.5331^, —0.5821, —0.6220, 
—0.7143, and —1.7610±0.1907z. That is, the approximate reduced quantum controller with four 
degrees of freedom also stabilizes the closed-loop system. 


4 Model Reduction for Completely Passive Linear Quantum 
Stochastic Systems 


Completely passive linear quantum stochastic systems are those which can be physically im¬ 
plemented using only passive components (that do not require external sources of quanta or 
energy); for example in quantum optics, these components are optical cavities, beam splitters. 


and phase shifters. In 44 , an asymptotically stable linear quantum stochastic system, in the 
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quadrature form 0 > is shown to be completely passive if and only if its controllability Gramian 
P is identity, i.e., P = I. The proposed choices of projection matrices Wq and Vq in the previous 
section does not guarantee that the controllability Gramian Pr of the reduced system in the 
quadrature form Q is P,. = / (when the reduced system is asymptotically stable). For this 
reason, we will now present a model reduction method for completely passive linear quantum 
stochastic systems, which guarantees both passivity and physical realizability properties of the 
reduced system, by appealing to the system dynamics described by annihilation operators. 

We begin by noting that each pair of position and momentum operators {qj,Pj) can be 
associated with a mode (or annihilation operator) as aj = {qj + ipj)l2, which serves as a 
quantum analogue of the amplitude of a harmonic oscillator. Let a = (oi, 02 ,..., a„)^ denote 
a vector of annihilation operators of a n degree of freedom quantum harmonic oscillators. The 
dynamics of a linear quantum stochastic system can be completely represented in the following 


annihilation-operator form if and only if the system is completely passive 14,16 


da{t) = Fa{t)dt + GdA{t); o(0) = a 
dY{t) = Ha{t)dt + KdA{t), 


( 20 ) 


where F G G G H G and D G Here, Y{t) and A{t) are as defined in 

Section [2j 

We will now present the physical realizability conditions of a linear quantum stochastic 
system in annihilation-operator form ( |20[ ). Similar to the earlier works [5,27 , when the number 
of outputs is less than inputs {i < m), we say that a completely passive linear quantum stochastic 
system is physically realizable if and only if there exists H G and K G such 

that the following augmented system (with I = m) is physically realizable: 

da{t) = Fa{t)dt -|- GdA{t)] a(0) = a 
H K 

dYait) = ^ a{t)dt + ^ dAit) (21) 

H \ [k 

where Ya{t) is a m-dimensional output vector. 

Theorem 9. A linear quantum stochastic system in annihilation-operator form ( |20[ ) (with 
i < m) is physically realizable if and only if 

F + F^ + = 0 ( 22 ) 

+ GK^ = 0 (23) 

and there exists K G such that the matrix K = (it'"'', iF is a complex unitary 

matrix in the sense that 

kk^ = k^k = L (24) 


Note that when i = m, K = K. 
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Proof. When £ = m, the necessity and sufficiency of (22)-(24) are shown in 14,16]. When 


i < m, the necessity and sufficiency of (22) and (24) follows immediately from the corresponding 
physical realizability conditions of the augmented system with £ = m. Let H = [H, H)~^ for 


some G £)xn^ For necessity of (23), we note from the corresponding physical realizability 


condition of the augmented system that = —GK"'. Post-multiplying both sides of this 


equation by , we have that H'^ = —GK^. For sufficiency of (23), let H = —KG^ 


where K is the complex matrix that makes (24) holds. From this H and (23), we have that 


= —GK^. This establishes the theorem statement. 


4.1 Reduced-order completely passive linear quantum stochastic systems 

Similar to the previous section, we will construct a reduced-order model of a completely linear 
quantum stochastic system via Galerkin projection. That is, given a full-order linear quantum 


stochastic system in annihilation form (20), we seek a reduced-order linear quantum stochastic 
system of the form 

dar{t) = Frar{t)dt + GrdA{t); 0^(0) = V^a 

dYr{t) = Hrar{t)dt + KrdA{t) (25) 

where F). = VaFVa, Gr = VaG, Hr = HVa, and Kr = K. Here, the projection matrix 14 G 
is orthonormal in the sense that 14^14 = F- 

We now propose our tangential interpolatory projection model reduction method. We will 
only present the left tangential interpolatory projection in this section, but the right tangential 
interpolatory projection can be similarly constructed involving a different subspace as demon¬ 
strated previously. To achieve the left tangential interpolation condition ([^ , consider a subspace 


Va = span I - F) ^ 


(26) 


I ' ' J i=l,2,...,r 

where interpolation points ui, ct 2 , ..., and the left tangent directions ni, ^ 2 , ■ ■ ■, fJ-r are chosen 
such that: (i) (cjjl — F) is invertible for all i = 1, 2,..., r, (ii) Hi A ^ alH = 1,2,..., r, and 
(iii) the subspace Va has the dimension r. Again, the condition on the dimension of Va ensures 
that the reduced-order system has r degree of freedom. The interpolation points and tangent 


directions can be chosen as discussed in Section 3.3, except that the points and directions do 


not need to be in conjugate pairs as I 4 can be a complex matrix. The following result is 
straightforward to show, so we shall simply state it without proof. 

Theorem 10. Given a full-order linear quantum stochastic system written in the annihilation- 


operator form (20) with the transfer function H(s), consider a reduced-order model of the form 


(25) with the projection matrix Va which is an orthonormal basis of the subspace Va defined 
in (26). Then the reduced-order model is physically realizable and its transfer function H,.(s) 


interpolates H(s) in the sense of ([^. 
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The tangential interpolatory projection method provides an avenue to reduce the order of 
a linear quantum stochastic system whilst ensuring that the reduced system is both physical 


realizable (as shown in Theorem 10) and completely passive (because the system dynamics are 
in annihilation-operator form). A key advantage of our method is that it is applicable to a 
larger class of completely passive linear quantum stochastic systems in comparison to some 
existing methods such as the quasi-balanced truncation method (which is applicable to 
those with unequal Hankel singular value of the product of the system’s controllability and 
observability Gramians, meaning that the system has less inputs than it does outputs), and the 
eigenvalue truncation method (which is applicable to systems with distinct eigenvalues). 
The following result is the analogue of Proposition for completely passive systems and is 
proved in an analogous way. 


Proposition 11. Consider a reduced-order model of the form (25) with the projection matrix Vo 


which is an orthonormal basis of the subspace Va defined in (26). LetUi{s) = ker{l/a (s/— F)^'} 


and U 2 {s) = ker{V'J(s/ — F)}. For any s G C that is not an eigenvalue of either F or F^, we 
have that 


|H(s) - H,(s)|| = \\H{sI - F)-Hi - Q(s))G|| , 
= ||F(/-F(s))(s/-F)-iG||, 


(27) 

(28) 


where Q{s) = {si — F)Va{sI — Fr)~^Va and R{s) = Va{sl — Fr)~^Va {si — F) are oblique 
projection operators (i.e., Q{s) and R{s) are idempotent), and when F and Fr have no poles on 
the right half plane and imaginary axis, we have the following Roo bounds: 

,-1/2 


^rlloo < sup 
ojeM 


1 - \\Pv± - P, 


Ui{iuj) I 


H{iujl - F)-^P^± 


P 




)G||) . 


< sup 
ojeM 


^ Pu2(^0J)\ 


- 1/2 


X HP 


1 X 2 ( 11 ^) I 


Pyx{zu;I - F)-^g\\) . 


(29) 


(30) 


4.2 Stability Property 

In this subsection, we will present some sufficient conditions that guarantee asymptotic stability 
of the reduced-order completely passive linear quantum stochastic system. 


Lemma 12. Given a linear quantum stochastic system written in annihilation-operator form 


(20), consider a reduced-order model of the form (25) with a projection matrix 14 that is an 


orthonormal basis of the subspace Va defined in (26). The reduced-order system {Fr, Gr,IIr, Kr) 
is asymptotically stable and is a minimal (controllable and observable) realization if and only if 


G^VaZr / 0 , 


(31) 
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for each eigenvector Zr of Fr. Moreover, (31) holds whenever one of the following sufficient 
conditions is satisfied: 


(i) Va C ker{Gt}^. 

(a) the interpolation points ai / lUJr for all i = 1,2 ,... ,r where lUJr is an eigenvalue of the 
resulting reduced-order system matrix F^. 


Proof. First note from Theorem (10) that the reduced system satisfies 

Fr + F^ + GrGl = 0 

Pre- and post-multiplying the above equation by zl and Zr, respectively, we have that 




Glzr 


= 0 , 


(32) 


where is the eigenvalue of Fr corresponding to Zr- This implies that the system is stable (i.e., 

= IlG^VaZrll / 0 for all eigenvectors Zr- 


5R{Ar} < 0 for all eigenvalues A^) if and only if 
Minimality of the reduced model then follows from 16, Lemma 2], where it is shown that for 


GrZr 


completely passive linear quantum systems asymptotic stability is equivalent to minimality. 
We now show the sufficiency of Conditions (i) and (ii). 

Sufficiency of Condition (i): Since 14 is an orthonormal basis of Va, we have that 14^;,. G Va- 
Now because Va 4 kerjil^}"'', G^VaZr / 0. 


Sufficiency of Condition (ii): First note from (32) that the poles of the reduced system do 
not lie on the right-half plane. We then prove the sufficiency of this condition by showing that 
when Ctzr = 0, the interpolation point is cTj = lUJr for some i = 1,2,... ,r. Note that this proof 


follows similar construction to the proof of 38, Theorem 11]. Suppose that Hj.(s) has a pole 
on the imaginary axis lUJr. Let Zr be the eigenvector corresponding to the eigenvalue lOJr of Fr, 
i.e., FrZr = uOrZr. Also note, from the physical realizability of the reduce system (established in 
Theorem 10), that Hr = —KrCt. Therefore, from this identity and (32), we have that HrZr = 0 


whenever Clzr = 0. 

Now let Vi = We also let V = [vi V 2 ... Vr] and U = [fii yi 2 ... ^rj- 

From the definitions of Vi, V, and U, we have that 


= 0 


(33) 


where S = diag((Ti, £72 ,..., Cr) is a diagonal matrix with the interpolation points on its diagonal. 
Note that V is also a basis of Va. Since the projection matrix 14 used to construct the reduced 
system is an orthonormal basis of Va, there exists a non-singular (full-rank) matrix Ta G 
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such that V = VaTa- Substituting VaTa for V and post-multiplying (33) by VaZr, we have that 


T^aV^FVaZr - ^T^VaZr + HVaZr = 0 
T^FrZr — '^T^Zr + Hj-Zr = 0 
lUJrTlzr — '^T^Zr = 0 

{lOJrlr — '^)Tlzr = Q (34) 


The 3rd step follows from the previously obtained result that H^Zr = 0. Since Ta has full rank 
and Zr is non-trivial (because it is an eigenvector), we have that T^Zr / 0. Moreover, because 
S = diag((Ti, (T 2 , ..., Ur), the result ( pd] ) implies that Uj = lUJr for some i = 1, 2,..., r. This 
establishes the sufficiency of Condition (ii) and hence the lemma statement. ■ 


We stress that the Condition (ii) of Lemma 12 is satisfied in many situations because it is 
generally improbable to find a purely imaginary interpolation point such that the point matches 
the imaginary pole of the resulting reduced system by chance. Hence, our proposed model 
reduction method for completely passive systems will generically lead to an asymptotically 
stable reduced-order linear quantum stochastic system. 


4.3 Illustrative example (passive system) 

Example 13. Consider a cascade connection of five identical optical cavities with the decay 
rate 7 = 10 ® and the resonant frequency ojq (uiq is much larger than 7 hut the exact value is 
not important here). Each cavity consists of two mirrors labeled Ml and M2. The mirrors Ml 
and M2 of cavity 1 are driven by coherent light fields Ai and M 2 , respectively. The amplitudes 
of these fields are modulated by a carrier frequency that is matched to ojq. The light reflected 
off Ml and M2 of cavity j drives the mirrors Ml and M2 of cavity j + 1, respectively. The 
two outputs of the system are the light reflected off the mirrors of cavity 5. Let Oj describes the 
oscillator mode of cavity j. Working in the rotating frame of reference with respect to wq, the 


dynamics of this system can be described in annihilation-operator form ( 20 ) with 


-7 

0 

0 

0 

0 


’-V 7 

-V7’ 

-27 

-7 

0 

0 

0 


-V 7 

-V 7 

-27 

-27 

-7 

0 

0 

, G = 

-V 7 

-V 7 

-27 

-27 

-27 

-7 

0 


-V 7 

-V 7 

.-27 

-27 

-27 

-27 



.-V 7 

-V 7 . 


H = —Gfi and K = F The frequency response (Bode plot) of the original system is illustrated 
in Fig. by black solid lines. Note that the response across the negative frequencies is a mirror 
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image of the response across the positive frequencies (i.e., the responses are symmetric around 
the origin) because F,G,H,K are real matrices. 

In this example, we are interested in obtaining an approximating system with r = 3. Note 


that the balanced truncation method in )21l is not suitable here because the product of the control¬ 


lability and observability Gramians is identity, leading to equal Hankel singular values of 1. The 
eigenvalue truncation method in l20j (truncating subsystems corresponding to larger eigenval¬ 
ues) is also not suitable because the poles of the system are all at —10®. Therefore, we apply the 
left tangential interpolatory projection model reduction approach. As suggested in Section 
we choose the tangent directions to be {p-i, fi 2 , tJ-s) = (ei,ei,ei). Note from the system dynamics 
that choosing 1 x 2 , p-s) = (ei,ei,ei) or {pi, p2, Ts) = {^ 2 ,^ 2 , £ 2 ) would result in the same 
subspace Va. Since the input-output responses of the original system have characteristics of both 
low-pass and high-pass filters, we simplify the optimization by assuming that the interpolation 
points are of the form (cxi, cr 2 , cjs) = {uv^,0, —ico^) (i.e., the interpolation point at 0 matches 
the low-frequency responses while the other points match the high-frequency responses). This 
form of interpolation points also preserve the symmetric properties of the frequency responses 
around the origin (i.e., F^, Gr, and Hr are real matrices). We find a local minimizer of the TL 2 
optimization problem to obtain the value of because the Hoc error is 2 for any > 0. We 
have that = 1.48 x lO'^. This leads to a stable system with the poles at —5.1541 x 10® and 
(-1.0780 ±0.8142z) x 10^. 

Fig. illustrates the frequency response (Bode plot) of the reduced system using red dashed 
lines. From the figure, we see that the frequency responses of our approximation are closely 
matched to those of the original system along the pass bands. The error ||H(s) — Hr(s)|| for 
each s is shown in Fig. and Again, the error is large when s is close to the poles of the 
systems. As previously mentioned, the Hoc error is 2.00. The Hoo error bounds (29) and (30) 


are both 2.92. Here, the bounds are conservative because passive quantum systems have bounded 
real transfer functions, i.e., ||H(za;)|| < 1 for all ui ^ M. which implies that the Hoc error 


- <2 
--r II00 — 


5 Conclusion 

In this paper, we have proposed a tangential interpolatory projection model reduction method 
for linear quantum stochastic systems. The proposed approach retains the required physical re¬ 
alizability property of the original full-order system while ensuring that the transfer function of 
the reduced-order system matches that of the original system at multiple points (or frequencies) 
along some tangent directions. That is, the reduced system can be designed to match specified 
input-output responses of the original system at various frequencies. We also establish an Hoc 
error bound for the proposed method, formulate optimization based routines to select interpola¬ 
tion points along the imaginary axis, and introduce a heuristic for selecting tangent directions. 
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Figure 6: Example 13 The frequency responses (Bode plots) of the original system (black solid line) 


and the tangential interpolatory projection approximation (red dashed line). 
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Figure 7: Example [T^ The approximation error ||S(s) — .^.^(s) 


A passivity preserving model reduction method has also been proposed. We establish a new 
result illustrating that our reduced-order completely passive linear quantum stochastic system 
will typically be asymptotically stable. Significantly, our tangential interpolatory projection 
approach can be applied to a wider class of linear quantum stochastic systems compared to 
existing model reduction methods for linear quantum systems. Several examples are provided 
to illustrate the merits of our proposed model reduction approaches. 
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Figure 8: Example 13 The approximation error for each purely imaginary point (or frequency) ||5(za;) — 
S^(za;)||. 
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